; docformat = 'rst'
;
;+
;
; :Purpose:
;   Combine Eulerian and Lagrangian model sensitivities
;   
; :Inputs:
;   sim_name: Simulation name
;   
; :Requires:
;   LPDM footprints in untis of (mol/mol)/m2/s::
;     {input_directory}/fp/{Station_code}_{year}.sav
;
;   EUM sensitivities::
;     {input_directory}/{sim_name}/EUM/sensitivity.inp.sav
;   
; :Outputs:
;   Creates {input_directory}/{sim_name}/sensitivity.sav file that contains 
;   sensitivity matrix and prior emissions
;
; :History:
; 	Written by: Matt Rigby, MIT, Aug 19, 2011
;
;-
pro cels_sensitivity, sim_name

  compile_opt idl2
  on_error, 2

  print, "CELS_SENSITIVITY: combining sensitivities"
  
  StartY_IC=cels_get_parameter(sim_name, 'STARTY_IC')
  StartY=cels_get_parameter(sim_name, 'STARTY')
  EndY=cels_get_parameter(sim_name, 'ENDY')
  Stations=cels_get_parameter(sim_name, 'STATIONS')
  Stations_all=[stations, cels_get_parameter(sim_name, 'BG_STATIONS')]
  molmass=cels_get_parameter(sim_name, 'MOLMASS')
  measure_Average=cels_get_parameter(sim_name, 'MEASURE_AVERAGE')

  restore, cels_filestr(/Input, sim_name + '/EUM/sensitivity.inp.sav')
  restore, cels_filestr(/Input, sim_name + '/state.sav')
  restore, cels_filestr(/Input, sim_name + '/regions.sav')

  ;Get spinup if StartY_IC is specified
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  
  if startY_IC ne !null then begin
    restore, cels_filestr(/Input, sim_name + '/measurements/measurements_spinup.sav')
    if total(MF_STATION_ID ne STATIONS_ALL) gt 0. then begin
      message, "EUM spinup stations don't match"
    endif
    mf_spinup=temporary(mf)
    mf_spinup_error=temporary(mf_error)
    mf_spinup_time=temporary(mf_time)
    mf_spinup_stations=temporary(mf_station_ID)
    mf_spinup_wh=temporary(mf_wh)
    mf_spinup_ai=temporary(mf_ai)
    restore, cels_filestr(/Input, sim_name + '/EUM/sensitivity_spinup.inp.sav')
    if total(EUM_SPINUP_STATIONS ne STATIONS_ALL) gt 0. then begin
      message, "EUM spinup stations don't match"
    endif
    if total(EUM_SPINUP_TIME ne MF_SPINUP_TIME) gt 0. then begin
      message, "EUM spinup time doesn't match measurements"
    endif
    nMeasure_spinup=n_elements(mf_spinup)
  endif else begin
    nMeasure_spinup=0L
  endelse

  ;Get measurements
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  
  restore, cels_filestr(/Input, sim_name + '/measurements/measurements.sav')

  ;Check that measurements are all aligned.  This could be removed if some proper checking was done instead.
  if total(EUM_STATIONS ne STATIONS_ALL) gt 0. then begin
    message, "EUM stations don't match"
  endif
  if total(MF_STATION_ID ne STATIONS_ALL) gt 0. then begin
    message, "Measurement stations don't match"
  endif

  nStations=n_elements(stations)
  nStations_all=n_elements(stations_all)

  ;Get regional total emissions
  restore, cels_filestr(/Input, sim_name + '/emissions/emissions_totals.inp.sav')
  if total(EUM_ID_IC ne pulse_magnitude_IC_ID) gt 0. then begin
    message, 'IC pulse magnitudes mixed up?'
  endif
  if total(EUM_ID_NLR ne pulse_magnitude_NLR_ID) gt 0. then begin
    message, 'NLR pulse magnitudes mixed up?'
  endif
  if total(EUM_ID_LR ne pulse_magnitude_LR_ID) gt 0. then begin
    message, 'LR pulse magnitudes mixed up?'
  endif


  ;EUM sensitivities
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

  ;Count aggregated regions
  N_AR=0L
  for lri=0, n_LR-1 do begin
    restore, cels_filestr(/Input, sim_name + '/ar_' + string(lri, format='(I02)') + '.sav')
    n_AR+=n_elements(ar_wh->keys())
  Endfor
  n_ART=n_AR*(EndY - StartY)  ;Annual only

  statesize=N_IC + N_NLRT + N_ART

  nMeasure=n_elements(mf) + nMeasure_spinup
  
  H=fltarr(stateSize, nMeasure)
  x_ap=fltarr(statesize)
  x_lr=intarr(stateSize)-999
  x_time=dblarr(stateSize)-999.d
  x_wh=hash()

  y=fltarr(nMeasure)

  ;IC sensitivity
  for xi=0, n_IC-1 do begin
    hTemp=reform(h_spinup[xi,*, *])
    H[xi, 0:nMeasure_spinup-1]=Htemp[mf_spinup_wh]
    hTemp=reform(h_IC[xi, *,*])
    H[xi, nMeasure_spinup:nMeasure-1]=Htemp[mf_wh]
  endfor
  x_ap[0:n_IC-1]=emissions_IC ;Not necessarily emissions, could be mole fractions MLR
  undefine, H_IC

  ;NLR sensitivity
  for xi=0, n_NLRT-1 do begin
    hTemp=reform(h_NLR[xi, *, *])
    H[n_IC+xi, nMeasure_spinup:nMeasure-1]=hTemp[mf_wh]
    year=fix(strmid(EUM_ID_NLR[xi], 2, 2))
    if year gt 50 then year+=1990 else year+=2000
    x_time[n_IC + xi]=julday(7, 1, year, 0)
  endfor
  x_ap[n_IC:N_IC+N_NLRT-1]=emissions_NLR ;kg/s
  undefine, h_NLR

  ;MEASUREMENTS
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  
  y=[mf_spinup, mf]
  y_station=[reform(mf_spinup_ai[0, *]), reform(mf_ai[0, *])]
  y_ai=[[reform(mf_spinup_ai)], [reform(mf_ai)]]
  y_time=[reform(mf_spinup_time[mf_spinup_ai[1, *]]), reform(mf_time[mf_ai[1, *]])]
  y_error=[mf_spinup_error, mf_error]
  caldat, y_time, y_month, y_day, y_year


  ;Coupled sensitivities
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  
  state_cumulative=n_IC + N_NLRT
  for yi=StartY, EndY-1 do begin  ;Annual only: no separation between y and x annual divisions
  
    if yi ge 2000 then Year_Str=string(yi - 2000, format='(I02)') else year_str=string(yi - 1900, format='(I02)')

    for lri=0, N_LR-1 do begin
      ;Get local emissions
      restore, cels_filestr(/Input, sim_name + '/emissions/lr_' + string(lri, format='(I02)') + '_' + string(yi) + '.sav')
      ;Get aggregated regions
      restore, cels_filestr(/Input, sim_name + '/ar_' + string(lri, format='(I02)') + '.sav')
      n_AR_LR=n_elements(ar_wh->keys())
      area=areagrid(lonq, latq)
  
      ;Aggregated emissions
      H_LR_index=where(EUM_ID_LR eq strcompress(string(lri, format='(I02)') + year_str, /remove_all), count)
      if count eq 0 then begin
        message, "Can't find correct EUM ID for background LR change"
      endif      
      H_LR_temp=reform(H_LR_BG[H_LR_index, *, *])
      for ari=0, n_AR_LR-1 do begin
        ar_indices=ar_wh[ari]
        x_ap[state_cumulative + ari]=total(q[ar_indices]*area[ar_indices])  ;kg/s
        
        ;Add background sensitivity
        H[state_cumulative + ari, nMeasure_spinup:-1]=H_LR_temp[mf_wh]
      
        ;State indentifiers
        x_time[state_cumulative + ari]=julday(7, 1, yi, 0)
        x_lr[state_cumulative + ari]=lri
        x_wh[state_cumulative + ari]=ar_wh[ari]
      endfor

      ;Combine footprints from multiple stations
      which_stations=where(station_lr eq lri)
      for si=0, n_elements(which_stations)-1 do begin
  
        ;locations of this year and station in y vector
        wh_time_y=where(y_station eq which_stations[si] and y_time ge julday(1,1, yi, 0) and y_time lt julday(1, 1, yi+1, 0), count)
              
        if count ne 0 then begin
        
          ;Get footprint
          restore, cels_filestr(/Input, 'fp/' + stations[which_stations[si]] + '_' + string(yi) + '.sav')
          
          ;Time indices where measurements are located in footprint file
          ;TIME STEPS IN FOOTPRINT FILE MUST BE THE SAME AS IN MEASUREMENT FILE
          wh_time_fp=mf_ai[1, where(mf_ai[0, *] eq which_stations[si])]
          dummy=min(abs(mf_time - julday(1, 1, yi, 1)), SYearI)  ;1 hr past midnight because first time step is 1.30am in NAME. MLR
          dummy=min(abs(mf_time - julday(12, 31, yi, 23, 59)), EYearI)
          wh_time_fp=wh_time_fp[where(wh_time_fp ge SYearI and wh_time_fp le EYearI)] - SYearI
          if n_elements(wh_time_fp) ne count then begin
            message, 'Something has gone wrong adding in LPDM footprints'
          endif
    
          ;Add aggregated footprints to background
          for ari=0, n_AR_LR-1 do begin
            ar_indices=ar_wh[ari]
            for ti=0, n_elements(wh_time_y)-1 do begin
              fp_temp=fp[*, *, wh_time_fp[ti]]/area/molmass*1000.*1.e12
              ;Aggregate footprints within each region as emissions-weighted average
              H[state_cumulative + ari, wh_time_y[ti]]+=$
                total(fp_temp[ar_indices]*q[ar_indices]*area[ar_indices])/total(q[ar_indices]*area[ar_indices])
            endfor
          endfor

        endif else begin
          ;Skip station, because there is no data in this year
        endelse
       
      endfor  ;which_stations
      state_cumulative+=n_ar_lr
  
    endfor  ;lr
  endfor  ;year

  ;Time-averaging
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  
  TimeSize=fix((julday(1, 1, EndY, 0) - julday(1, 1, starty, 0))/double(Measure_average))
  time=julday(1, 1, StartY, 0) + dindgen(TimeSize)*double(Measure_average) + double(Measure_average)/2.

  y_ti=lonarr(nMeasure - nMeasure_spinup)
  
  For ti=0, nMeasure - nMeasure_spinup-1 do y_ti[ti]=fix((y_time[nMeasure_spinup + ti]-julday(1, 1, StartY, 0))/double(Measure_average))
  wh=where(y_ti ge timeSize-1)
  if wh[0] ne -1 then y_ti[wh]=timeSize-1

  sortstr=strcompress(string(y_ti, format='(I04)') + string(y_station[nMeasure_spinup:-1], format='(I02)'), /remove_all)
  sortstr=sortstr[uniq(sortstr, sort(sortstr))]
  nMeasure_new=n_elements(sortstr)
  y_tiu=long(strmid(sortstr, 0, 4))
  y_ti=[replicate(-999, nMeasure_spinup), y_ti]

  y_time_new=dblarr(nMeasure_spinup + nMeasure_new)
  y_error_new=dblarr(nMeasure_spinup + nMeasure_new)
  y_new=fltarr(nMeasure_spinup + nMeasure_new)
  h_new=fltarr(StateSize, nMeasure_spinup + nMeasure_new)

  if nMeasure_spinup gt 0 then begin
    y_station_new=[y_station[0:nMeasure_spinup-1], long(strmid(sortstr, 4, 2))]
    y_time_new[0:nMeasure_spinup-1]=y_time[0:nMeasure_spinup-1]
    y_new[0:nMeasure_spinup-1]=y[0:nMeasure_spinup-1]
    y_error_new[0:nMeasure_spinup-1]=y_error[0:nMeasure_spinup-1]
    y_station_new[0:nMeasure_spinup-1]=y_station[0:nMeasure_spinup-1]
    h_new[*, 0:nMeasure_spinup-1]=h[*, 0:nMeasure_spinup-1]
  endif

  For n=0, nMeasure_new-1 do begin
    y_wh=where(y_station eq y_station_new[nMeasure_spinup + n] and y_ti eq y_tiu[n], count)
    if count eq 0 then stop
    y_new[nMeasure_spinup + n]=mean(y[y_wh])
    y_time_new[nMeasure_spinup + n]=time[y_tiu[n]]
    for xi=0, StateSize-1 do h_new[xi, nMeasure_spinup + n]=mean(h[xi, y_wh])

    ;Error calculation
    station_temp=y_station_new[nMeasure_spinup + n]
    if station_temp lt nStations then begin
      lr_temp=station_lr[station_temp]
      h_temp=reform(h[where(x_lr eq lr_temp), *])
      h_temp=h_temp[*, y_wh]
      x_temp=x_ap[where(x_lr eq lr_temp)]
      y_error_new[nMeasure_spinup + n]=sqrt((stddev(h_temp##x_temp, /nan)/sqrt(float(Measure_average)/5.))^2 + mean(y_error[y_wh])^2)
    endif else begin
      ;Background station
      y_error_new[nMeasure_spinup + n]=sqrt(mean(y_error[y_wh])^2)
    endelse
    
  endfor

  y_station=temporary(y_station_new)
  y_time=temporary(y_time_new)
  y_error=temporary(y_error_new)
  y=temporary(y_new)

  h=temporary(H_new)
  
  save, filename=cels_filestr(/Input, sim_name + '/sensitivity.sav'), $
    h, x_lr, x_time, x_ap, y, y_time, y_station, y_error
  
End